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I. INTRODUCTION 

Complex network theory represents a general unify- 
ing formalism under which is possible to understand and 
rationalize the intricate connectivity and interaction pat- 
terns of many natural and man-made systems. Thus, a 
systematic statistical analysis of large scale dataset has 
allowed to unveil the existence of apparently universal 
topological features, shared by a large number of differ- 
ent systems from the technological, social or biological 
domains [TU5]. Among these characteristics, probably the 
most intriguing is the discovery of the apparently ubiqui- 
tous scale-free (SF) nature of the connectivity described 
by a probability P{k) ~ that an element (vertex) 
is connected to other k elements (has degree k), with a 
degree exponent usually in the range 2 < 7 < 3 [U [3] . 
These and other discoveries have promoted a large mod- 
eling activity, aimed at understanding the origin and na- 
ture of the observed topological features |2j [3| . In recent 
years, the research community has also devoted a great 
deal of attention to the study of the dynamical processes 
on complex networks |?l ISl, which can have important 
implications in the understanding of real processes such 
as the spread of epidemics in social systems [6] or traffic 
in technological systems as the Internet [71 or transport 
infrastructures [8^. 

The theoretical understanding of dynamical processes 
on complex networks is based in the application of mean- 
field approaches that are essentially based in the annealed 
network approximation [S]. Any non- weighted network 
is fully characterized by its adjacency matrix 0^-, taking 
the value 1 when vertices i and j are connected by an 
edge, and zero otherwise. In real (quenched) networks, 
the values of the adjacency matrix are fixed and do not 
change with time. When a dynamical process takes place 
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on top of such a network, one is considering the network 
as frozen, with respect to the characteristic time scale tjj 
of the dynamics. Other networks, however, are dynam- 
ical objects, changing in time over a characteristic time 
scale Tat. In this case, the adjacency matrix is defined 
only in a statistical sense, and a complete description of 
the network can be given in terms of its probability dis- 
tribution P{k) and its degree correlations P{k\k') [9.. In 
the limit <ti td, when all edges are completely reshuf- 
fled between any two dynamical time steps, the resulting 
network is called annealed [101. Annealed networks rep- 
resent an extremely important theoretical tool, because 
mean-field predictions for dynamical processes turn out 
to be exact in this kind of substrates [51. 

Dynamical processes with absorbing configurations 
constitute a subject of outstanding interest in non- 
equilibrium Statistical Physics [HI [T2] that have also 
found a place in network science, as representative mod- 
els of practical problems ranging from epidemic spreading 
[13l[T4], infrastructure's resilience [El [16], or activated 
dynamics [17] . to mention just a few. The simplest lattice 
model allowing absorbing configurations is the classical 
contact process (CP) [H], whose universality class and 
mean-field (MF) description have been discussed along 
the last few years [TUl fl9 - 22j . In the CP defined in an 
arbitrary network, vertices can be in two different states, 
either empty or occupied. The dynamics includes the 
spontaneous annihilation of occupied vertices, which be- 
come empty at unitary rate, and the self-catalytic oc- 
cupation of an empty vertex i with rate Xrii/ki, where 
rii is the number of occupied neighbors of i and ki is 
its degree. The model is thus characterized by a phase 
transition at a value of the control parameter A = Ac, 
separating an active from an absorbing phase devoid of 
active vertices. Despite of its simplicity, the CP on SF 
networks exhibits a very complex critical behavior even 
if studied in an annealed substrate [101 [HI [12] . This case 
is particularly interesting, since in annealed networks all 
connections are rewired at a rate much larger than the 
typical rates involved in the dynamical process, implying 
that dynamical correlations are absent ITU!. In this case. 
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the MF approach is expected to be an exact description 
of the problem. 

The configuration in which all vertices are empty plays 
a very particular role, since once the system has fallen 
into this state, the dynamics becomes frozen. For this 
reason, these states are called absorbing and constitute 
a central feature in the analysis of finite size systems 
since, in this case, the single actual stationary state 
is the absorbing one [TT]. Finite size and absorbing 
states must therefore be handled using suitable strate- 
gies, concomitantly with an ansatz for the finite size scal- 
ing (FSS) [23]. A widely adopted procedure is the so 
called quasi-stationary (QS) state [MIHSIj in which the 
absorbing configuration is suitably excluded from the dy- 
namics. 

In this work, we present a study of the QS state of the 
CP on SF annealed networks, combining the QS numer- 
ical approach developed in Ref. suitably extended 
to complex networks, with the theoretical analysis of an 
approximated one-step process derived from mean-field 
theory [2T] . Our analysis allows us to determine the prob- 
ability distribution of activity both close to the critical 
point and in the off-critical regime, as well as to obtain 
high quality data for relevant QS quantities, such as the 
density of active sites or the characteristic times. This 
last information is used to check the finite-size scaling 
forms derived from mean-field theory, which turn out to 
be loaded with very strong corrections to scaling. 

The results presented in our paper provide a deeper 
understanding of the nontrivial dynamics of the contact 
process in annealed networks. Moreover, they open the 
path to extension of the QS approach to the analysis of 
other dynamic processes with absorbing phase transitions 
taking place on more complex and/or realistic substrates 
such as, for example, quenched networks and small- world 
topologies |26| . where edges are never rewired and dy- 
namic correlations are usually present. 

We have organized our paper as follows. In Sec. [ll] 
we review the necessary background for QS analysis and 
numerical simulations, while in Sec. |III| we summarize 
the main results of MF theory for the CP on annealed 
networks. A master equation approach for the QS state 
is developed in Sec. |IV| Section |V] is devoted to discuss 
the finite-size scaling forms of the relevant QS quantities, 
as well as the corrections to scaling at criticality. The 
off-critical analysis is discussed in Sec. |VI[ F inally, our 
concluding remarks are presented in Sec. VII| 



II. FINITE SIZE AND THE 
QUASI-STATIONARY STATE 

In finite systems, the absorbing state is a fixed point 
that can be visited even in the supercritical phase due to 
stochastic fluctuations. Numerical simulations of finite 
systems are particularly sensitive to absorbing states and 
therefore suitable simulation strategies are required. The 
standard procedure consists in restricting the averages to 



those runs that did not visit the absorbing configuration 
[llj . the so-called surviving averages. From a mathemati- 
cal point of view, it is useful to define the quasi-stationary 
(QS) state that consists of the ensemble of states accessed 
by the original dynamical process at long times restricted 
to those not trapped into an absorbing one [24] . The in- 
tensive quantities in a QS ensemble must converge to the 
stationary ones in the thermodynamic limit. Thus, in 
the active phase, the lifespan grows exponentially fast 
with the system size and the QS state becomes identi- 
cal to the stationary one. In the sub-critical phase, on 
the other hand, the activity in the QS state corresponds 
only to a few [0(1)] particles fiuctuating above the ab- 
sorbing state implying a density that vanishes inversely 
proportional to the system size. 

Formally, the QS state is related with the original one 
in the limit t ^ cx) by 



(1) 



where P{cr) is the QS probability associated to the state a 
and Ps{t) is the survival probability, i.e., the probability 
that the system is active up to time t. For a one-step pro- 
cess, the state of the system is completely determined by 
the number of occupied vertices n. Letting Pn(t) be the 
probability that the system has n particles at time t, the 
QS distribution P„ is given by P„(t) — Ps{t)Pn, for which 
the normalization X]n>i Pn = ^ applies. The probability 
of visiting the absorbing state is redistributed among the 
active configurations, proportionally to P„, which consti- 
tutes the essence of the QS state [24]. Knowledge of the 
QS distribution allows to compute the standard quanti- 
ties associated to this state. For example, the probability 
to visit the vacuum in the CP is given by Po = Pi , inde- 
pendently of the network substrate. Thus, it is straight- 
forward to show that the survival probability and the 
pre-absorbing state are related by dPg/dt = — PiP, pro- 
viding a characteristic time scale 



1 

P~i' 



(2) 



Analogously, the density of active sites in the QS state is 
given by 



- y np„ 



(3) 



The standard numerical procedure to simulate the QS 
regime based on averages over survival runs has a limited 
accuracy due to the very rare achievement of surviving 
configurations at very large times. Thus, for instance, 
the stationary densities are determined as a plateau at 
long times in the curve p(t) that is usually noisy and 
short close to or below the criticality due to the limited 
number of independent runs computationally accessible. 
The previous interpretation of the QS state provides an 
alternative simulation strategy, in which every time the 
system visits the absorbing state, this configuration is re- 
placed by an active one randomly taken from the history 
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of the simulation [53] . For this task, a list with M active 
configurations is stored and constantly updated. An up- 
date consists in randomly choosing a configuration in the 
list and replacing it by the present active configuration 
with a probability pr- After a relaxation time tr, the QS 
distributions are determined during an averaging time ta- 
The improved QS method has been successfully applied 
to accurately determine the universality class of several 
models with absorbing configurations [211 dZl HH] • 



III. CP IN ANNEALED NETWORKS 

The network in which dynamics takes place is assumed 
to be annealed and therefore completely characterized by 
the degree distribution P{k), the probability that a ran- 
domly chosen vertex has a degree k, and the degree cor- 
relation function P{k'\k) defined as the conditional prob- 
ability that a vertex of degree k is connected to a vertex 
of degree k' [5]. The number of vertices of the network 
is denoted by N and its maximum degree (cutoff) by kc 
[29] . In an annealed framework, the MF rate equation for 
the density of occupied vertices in the degree class k (i.e. 
the probability that a vertex of degree k is occupied) is 
given by (21j 



d , , 



P{k'\k)pk' 



(4) 



Close to the criticality, when the density at long times 
is sufficiently small such that pkc ^ 1, equation ([t]) be- 
comes pk — \kp/{k). Substituting in the transition rates, 
the first order approximation for the one-step process is 



W{n — l,n) — n, 

W{n+l,n) = Xn{l-Xgn/N), 



(8) 



in which g = (fc^)/(fc)^. Based on numerical evidences 
and scaling arguments, later confirmed by more rigor- 
ous means in Refs. [TOl HI], the authors proposed that 
the critical characteristic time r and stationary density 
p scale as ^ 



1/2 



and 



'1/2 



(9) 



(10) 



respectively. For a network with degree exponent 7 and 
a cutoff scaling with the system size as kc ~ N^^'^ , where 
Lo is an arbitrary positive parameter, the factor g scales 
for large N as g ^ for 2 < 7 < 3 and g ^ const 

for 7 > 3. Therefore the critical QS density scales as 
p ^ N^" where 



. 1 /3-7 

V = - + max — — , 
2 \ 2uj 



(11) 



The first term represents the spontaneous annihilation 
and the second one the creation inside the compartment 
k due to the interaction with all compartments under the 
hypothesis that there are no dynamical correlations. A 
simple linear stability analysis allows to show the pres- 
ence of a phase transition, located at the value Xc — 1 and 
independent of the degree distribution and degree cor- 
relations, separating an active from an absorbing phase 
with pk = [lOj . Considering in addition uncorrelated 
networks with P{k'\k) = k'P{k')/{k) [2], the overall den- 
sity p = pkP{k) obeys the equation 



dt 



p{t) = -p{t) + Xp 



i-{kr^Y.^p{k)pk{t) 



(5) 



In Ref. [21] , it was realized that the low density regime 
ofEq. m can be understood as a one-step process (bi- 
ased random walk) with transition rates 



Similarly, the characteristic time follows t ^ N " with 
exponent 



«=2-maxf— ,0 



(12) 



The MF supercritical density for an infinite system was 
found to vanish at criticality as p ~ where (3 = 1/(7— 
2) and A = A — Ac- For finite systems, the QS density 
has an anomalous cutoff-dependent FSS given by [T^ 



p(A,iV) = 




for 



A X(Jt)_ 
g kc 



(13) 



where G{x) ^ x for a; 3> 1 and G{x) is constant for 
X <C 1. The anomaly lies on the supercritical density 
dependence on the system size through the factor g given 
by A/5 if A>yg7]V[IU]. 



W{n — 1, n) 
W{n^ l,n) 



An [l-{k)-^Y.kkP{k)pk{t)\ 



(6) 



where W{m, n) corresponds to the transition from a state 
with n particles to one with m particles. The stationary 
state dtPk =0 of Eq. Q reads as 



Pk 



Xkp/{k) 

1 + xkp/{k) ■ 



(7) 



IV. MASTER EQUATION APPROACH OF THE 
QS STATE 

In order to gain analytical information about the QS 
distribution not far away from the critical point, we can 
consider the one-step process approximation described 
by the transition rates in Eq. (|6|. Starting from them, 
it is possible to write down a master equation (ME) for 
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the evolution of the number of particles P„ (t) , taking the 
standard form 

Pn = J2 ^("' ^)Prnit) - ^ W{m, n)P,,{t). (14) 

m m 

In the long time limit we have « and, consequently, 
equation ^ can be applied resulting the ME 

P„ = (n + l)P„+i + u„_iP„_i - (n + u„)P„, (15) 
with Un ~ Xn{l — Q) and Q given by 



e[p] 



El 



k^P{k) 
+ \kp/{k) ■ 



(16) 



lows !25_ 
to find Pi°\ n 



Substituting now P„(t) = Ps{t)Pn and using dPg/dt = 
— PiPg [25], the following recurrence relation is obtained 

Pn = -[(Wn-l + n - 1 - Pl)Pn-l " U«-2Pn-2], (17) 

n 

where n = 2, ■ ■ ■ , N and Pq = 0. The QS distributions 
are compl etely determined since Pi, the initial condition 
to iterate (17), is given by the normalization J2n>i = 
1. 

Full information of the QS distribution can be ob- 
tained from the ME by solving it numerically. A nu- 
merical recipe to iterate the recurrence relation is as fol- 
Start with a guess for P^^°^ and iterate (17) 
2, • • • ,N. Repeat the procedure using 
pO'+i) _ P^-'^ I^^P^'' as new guess until the normal- 
ization be reached. Suitable truncations can be used to 
speed up the numerical process and to prevent instabil- 
ities. The truncations at finite densities are justifiable 
since the Central Limit Theorem guarantees that fiuctu- 
ations much larger or much smaller than the average are 
exponentially negligible. 

In order to explore the properties of the QS state, 
we have performed extensive Monte Carlo simulations 
of the CP on annealed networks. We use a random 
neighbor network (RNN) [10] with degree distribution 
P(fc) - k-1, degree correlations P{k'\k) = k'P{k')/{k), 
a degree cut-off kc — N^/'^ , and a fixed minimum de- 
gree fco = 2. The single effect in increasing the minimum 
degree ko is a shift to higher densities which does not 
affect the critical properties. In an annealed approach, 
all links are redefined between any two time steps in such 
a way that the neighbor of a given vertex is selected by 
randomly choosing a vertex of the network with a prob- 
ability k' P{k') / {k) . We perform stochastic simulations 
using the usual scheme jllj : At each time step, an occu- 
pied vertex is chosen at random and the time updated as 
t ^- t -f At, where M = 1/[(1 -I- A)n(i)] and n{t) is the 
number of occupied vertices at time t. With probabil- 
ity p = 1/(1 + A), the occupied vertex becomes vacant. 
With complementary probability 1 — p = + X), one 
of its neighbors (following RNN rules) is selected and, if 
empty, occupied. If the selected neighbor is already oc- 
cupied the simulation goes to the next step. QS states 
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FIG. 1: (Color on-line) QS probability distributions at criti- 
cality obtained in simulations of the CP on annealed SF net- 
works with a cutoff N^^'^ are compared with the numerical 
solutions of Eq. (171. In the main plot the networks have size 
N = 2 X 10*. Inset A shows the curves of the main plot in a 
logarithmic scale in order to compare the tails. Inset B shows 



the QS distributions for a smaller network of size A'^ 
with 7 — 2.75. 



10^ 



were simulated using the method described in Sec. |TT] 
with M = 10^, pr = 0.02 At and ta = U = 10^. Due the 
short distance between vertices, the relaxation times are 
very short if compared with critical relaxation on regu- 
lar lattices. The network sizes were varied from 10'^ to 
10^ and 50 — 500 network samples were used in the aver- 
ages (the larger the size of the network, the smaller the 
number of samples). During the averaging interval, the 
current configuration is counted in the QS distribution 
with a probability proportional to its lifespan. 

The QS probability distributions for obtained in simu- 
lations of the critical CP are shown in Fig. [ijfor different 
values of the degree exponent 7 and network size N and 
compared with the results of the numerical solution of the 
recursion relation, equation (17 1. A remarkable agree- 



ment between the simulations and the one-step process 
approach is achieved even for the asymptotic (Gaussian) 
tail as one can see in the inset A of Fig. [l] A good ac- 
cordance, which is improved as size increases, is observed 
even for sizes as small as 10'^, while neat discrepancies ap- 
pears for N ~ 10^ (inset B of Fig. [l]). Actually, we can 
show that the one-step and the Langevin approach devel- 
oped in Ref. jlOj are equivalent in the low density limit 
(see Appendix) and, consequently, the lower the density 
the better the one-step mapping. 



V. THE QS STATE AT CRITICALITY 

A. Analytical approximation at criticality 

At criticality, where the densities at long times are very 
low, we have Un = An(l — An/fi) where f2 = N/g- In this 



limit, equation (15 1 corresponds exactly to the ME of 



the CP on a complete graph of size 17, for which the QS 
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analysis was already worked out elsewhere [5S] . Analyti- 
cal insights about the criticality can be obtained through 
a van Kampen's expansion |30j of the recurrence relation 
(17). Let us consider the scaling solution of Eq. (171 



0.008 



1 



f 



(18) 



where f{x) is a scaling function to be determined. Plug- 
ging Eq. ( |18[ ) into (17 1, and performing a Taylor expan- 
sion up to second order, the result up to order 51^^ is 



+ {2 + x^)^+2xf ^-M (19) 



where /o = Pi^l^/^ — /(O) must be chosen to impose the 
normalization condition f{y)dy = 1. Dickman and 
Vidgal ^25j analyzed Eq. ( 19 ) numerically and checked 



the agreement with the recurrence relation for the CP 
on the complete graph. We complement the analysis 
by obtaining the asymptotic behaviors analytically. It 
is straightforward to see that the distribution decays lin- 
early as f{x) ~ /o(l — fox/2) for x <C 1. A correction to 
this initial behavior can be obtained discarding the term 
xf"{x) (a low curvature approximation) and the solution 
satisfying the boundary condition /(O) = /o is 



fix) 



2/o 



cxp 



A 

del 
2/ 



arctan 



a; < 1. (20) 



For X 3> 1 the zeroth order terms are discarded and Eq. 
p9| turns to /" + xf + 2f ~0. The solution satisfying 
the boundary condition xf{x) — )• for x — >■ cxi {{x) is 
finite) is 



/(x) ~exp(-a;V2), a; > 1, 



(21) 



implying a Gaussian tail. 

In Fig. [2] we compare the results of numerical simu- 
lations of the CP at criticality (A = 1) on a network 
with N — 160000 nodes, degree exponent 7 = 2.25 and 
cutoff scaling exponent u — 2 with the corresponding 
analytical approximations. Despite of the lack of rigor, 
the analytical result is in good agreement with numerical 
simulations. The approximation (20) agrees with numer- 
ical analysis even for a; w 1 or equivalently for a number 
of active vertices n « \/f2. For the particular network of 
Fig. [2] we have « 147. The accordance still holds for 
small networks (^ 10"^) independently of the exponent 
degree and cutoff scaling. 



The scaling function (18) encloses the FSS form of the 
CP at criticality. In fact, the mean number of occupied 
vertices at the QS regime is given by 



N 



N 



(22) 



where f2 = N/g. Letting x = n/\/Cl and Ax = 1/vil, 
the sum can be approximated by a continuous integration 
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FIG. 2: Comparison between the numerical simulations of 
the CP on annealed networks at criticality and the asymp- 
totic behaviors of the QS distribution given by Eqs. (201 and 
(2l]). Numerical simulations for iV = 160000, 7 = 2.25 and 
(jj = 2 are represented by circles and the asymptotic solu- 
tions by lines. The dashed line is the linear approximation 
for n/v^ <^ 1. The inset shows the comparison with the 
Caussian tail using a straight line of slope —1/2(1. 



when N ^ 00, namely 

n = J2 xf{x)Ax w / _ ^1/2 _ 

x—m/ 

(23) 

Therefore, the critical QS density is simply p = n/N ^ 
{gN)~^/^, recovering the result firstly presented in 
Ref. [21]. Analogously, the characteristic time r is also 
directly obtained by the present QS analysis by noticing 
that, since A = /of^^^/^, we have t = 1/Pi ^ fl^^^ = 
{N/gY/\ 



B. Scaling at criticality 

The analysis of the QS state, either by direct QS sim- 
ulations or of the iterative solution of the corresponding 
approximate ME, allows us to obtain high-quality data 
for the characteristic quantities at criticality, namely the 
stationary density p and the characteristic time r. At 
criticality, these quantities are expected to exhibit a scal- 
ing with system size of the form p ~ N~'^ and r ~ iV", 



with exponents given by Eqs. ( 11 ) and ( 12 ) and depend 



ing on the cutoff scaling exponent uj. 

QS critical densities as functions of the network size 
iV, computed for several degree exponents and lo = 2 are 
shown in the main plot of Fig.|3] A gain, an incontestable 
agreement between QS simulations and the numerical 
ME approach is observed. Similar agreement is obtained 
for the analysis of characteristic time (data not shown) . 
However, if one tries to recover the theoretical scaling 
exponents, taking the values v = max[l/2, (5 — 7)/4] and 
a = min[l/2, (7 — l)/4], by means of a direct power law 
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FIG. 3: QS critical density for the CP on annealed SF net- 
works with cutoff kc = N^^^. In main plot, the densities 
obtained with the numerical solution of the ME (lines) and 
QS simulations (symbols) are shown. Inset A shows the ex- 
ponent of the scaling law p ~ A''"" obtained analytically for 
asymptotically large systems {N — >■ oo), performing a direct 
fit to a power-law form (PL) of the data in main figure, and 
performing a fit to a power with corrections to scaling. Inset 
B shows the same analysis for the characteristic time r ~ A''" . 



regression (insets in Fig. |3| , a very poor agreement with 
the expected analytical exponents is observed, as already 
noted in Ref. |31j. Indeed, if we look carefully at the 
data we can observe that, even though a pretty good lin- 
ear fit can be resolved for data range corresponding to 
10^ < < 10'^ as shown in Fig. H) the actual regressions 
are a little bit curved. Indeed, a careful analysis of the 
numerical data can resolve a slight downward (negative) 
curvature for 7 < 5/2 and an slight upward (positive) 
curvature for 7 > 5/2 at A = Ac. In a QS analysis these 
behaviors usually indicate a system slightly out of the 
critical point, being sub and supercritical for down and 
upward curvatures, respectively. But this is not the case 
for the data shown in Figs. |3] and |4] since the MF result 
is exact for the critical CP on annealed networks. Similar 
behaviors occur for plots of r versus N. 

As noted in [lOl [21] , the origin of this poor agreement 
between theory and simulations lies in the implicit de- 
pendence on N of the g factor defining the size scaling of 
p and T. In fact, the scaling forms p ~ 7V~^ and r ~ TV" 
make only sense in the limit of very large TV, when g 
has achieved its truly asymptotic form. For intermediate 
values of N, instead, one should keep the scaling forms 
with the simultaneous dependence on g and N [21j . If we 
want instead to make explicit the scaling with network 
size, we must consider that g ~ (fc^)/(A;)^ behaves, in the 
continuous degree limit, as 



(7-2)X"'' {i-c-')a-e-') 

(7-l)(3-7) (l-e^-2)2 



kl-\ (24) 

^(3-7)/w 



where ^ — ko/kc < 1. So, when N — >■ 00, g 
for 2 < 7 < 3 and g ^ const for 7 > 3. 

From Eq. ([24| it is possible to work out the explicit 



FIG. 4: Power law regressions (solid lines) of the critical QS 
densities (symbols) for CP on annealed SF networks with cut- 
off fee = N^^^ and degree exponents 7 = 2.25 and 7 = 3.00. 
Inset shows the QS densities for 7 — 2.25 rescaled by pure 
power law (PL), — pN^'^^'^^ , and a PL with correction to 



scaling, C = piV"-''**^^(l + 2 X 2"-^-'A/- 



' 4-3x2"-=' AT - 




FIG. 5: Critical QS characteristic times for CP on annealed 
SF networks with cutoff = A^i/It-i), 

main plot, lines 

represent numerical solution of the ME and symbols QS sim- 
ulations. The inset A shows the exponent of the scaling law 
p ~ A''^'^ while the inset B shows the exponent of the scaling 
law r ~ A''". Legends as in Fig. [s] 



form of the corrections to scaling in a direct analysis of 
the QS quantities as functions of the network size. So, 
considering 2 < 7 < 3 and performing an expansion to 
leading order in equation (24) yields 



9 ' 



const X (1 — ^ 



3-7 



2e 



7-2 , 



u3-7 



(25) 



Substituting into the MF scaling result p ^ (gN) 
we obtain the expression for the stationary density 



lnp^C-i)\nN+l'%^ 
^ N -J 



(26) 



Notice that the corrections do not introduce any parame- 
ters to be fitted. Similar expressions are found for 7 > 3. 
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Equation ( 26 ) explains the deviations from the of power 
law regime p cx observed for the CP on annealed 

SF networks. It is easy to see that the leading term for 
2 < 7 < 5/2 is negative and causes a downward curva- 
ture, in the same way that the leading term for 7 > 5/2 
bends the curve upwardly. Even though the correction 
vanishes for — > 00, it may occur extremely slowly due 
to the small exponents involved. For 7^3 and 7 2 
the corrections are logarithmic and are thus relevant for 
any finite size. 

Introducing the corrections given in Eq. (261 in the 
form 



In p' = In /9 



1 K 



2n' 



C-iylnN (27) 



and performing a linear fit, the asymptotic exponent {> 
is recovered as one can see in the inset A in Fig. [3] 
Equivalent corrections can be easily obtained for Inr vs. 
In N and the expected exponent a is recovered as shown 
in the inset B in Fig. |3] Additional proof of the strong 
finite size corrections is provided in the inset of Fig. [4] in 
which the critical QS density is rescaled by the predicted 
power law (PL) with exponent 9 and by this same PL 
with the correction (26). The first case is clearly size 
dependent while the second is flat. It is worth to note 
that the corrections are so strong for 7 = 2.25 and lu = 2 
that keeping only the leading term C'(A~"-^^^) was not 
enough to account for the deviation. 

We additionally performed the analysis for a; = 7 — 1, 
which corresponds to the natural cutoff that emerges 
in the absence of a structural cutoff [S]. For sake of 
simplicity, a hard cutoff was adopted such that connec- 
tivities larger than ~ 7Vi/(7-i) ^re forbidden. For 
this cutoff, the scaling exponents for the critical density 
and characteristic time are v — max[l/2, 1/(7 — 1)] and 
a = max[l/2, (7 — 2)/(7— 1)], respectively. Exactly as in 
the w = 2 case, the QS analysis via ME agrees with sim- 
ulations and the correct scaling exponents are obtained if 
corrections to the scaling are considered. Figure [5] shows 
the characteristic time and the insets therein the scaling 
exponent analysis for w = 7 — 1. Comparing the insets 
m Figs. |3] and [5) we can observe that the relative devia- 
tion between the exponents obtained using a simple and 
a corrected PL is smaller for the natural (w = 7 — 1) than 
for the cutoff kc = N^^^. Indeed, equation (26) tells that 



the larger the cutoff exponent oj, the stronger the cor- 
rections to the scaling. However, even if the cutoffs are 
not imposed, the natural one emerges spontaneously in 
networks with power- law degree distributions [5]. Con- 
sequently, these corrections to the scaling may be also 
present in the CP and other dynamical processes in SF 
substrates including the quenched case. 

The previous results may have a remarkable impact in 
the analysis of absorbing phase transitions in complex 
networks. The usual QS analysis assumes a power law 
dependence of the order parameters with the size at crit- 
icality. Such assumption is commonly used as a criterion 




FIG. 6: (Color on-line) Supercritical densities as functions 
of the system size for three degree distributions and a cutoff 
kc = N^^^. Solid lines are the ME numerical solutions and 
symbols QS simulations. 



to determine the critical point of absorbing phase tran- 
sitions in regular lattices [TT] and has been extended to 
quenched complex networks [121 [50] . Corrections to scal- 
ing in the form 1 — const x iV~°-^^ were already observed 
in QS simulations of the directed percolation universality 
class in hypercubic lattices, including the contact process 
[32] ■ Since these corrections decay with a large ex- 
ponent, they are significant only for small systems. In 
complex networks, the scenario is quite different since the 
corrections, which emerge from the intrinsic SF nature of 
the substrate, vanish very slowly and are important even 
for large systems {N ^ 10^ in the present work). 



VI. OFF-CRITICAL QS ANALYSIS 

The analysis performed for the critical CP in the pre- 
vious sections is expected to work also for the off-critical 
phase if the densities are still sufficiently small. In Fig. 
|6]we compare the QS density obtained by ME iterative 
solution and numerical computed for networks with dif- 
ferent degree exponents. In the plot we explore the su- 
percritical regime with rates A — 1.004 and 1.040 (0.4% 
and 4% above the critical point, respectively). The one- 
step ME predicts QS densities very accurately even for a 
substantial distance from the critical point, corroborat- 
ing that the approach is also suitable for the supercriti- 
cal phase. The anomalous FSS form in the supercritical 



regime, equation (13), depending on A, iV and g simul- 
taneously, is checked in Fig. [7j where we present the 
collapses of the data obtained from iterative solutions of 
the ME. Degree exponents 7 = 2.25 and 2.75 using struc- 
tural (a; — 2) and natural (w = 7 — 1) cutoffs are shown. 
Excellent collapses arc obtained in all cases, in agreement 
with Refs. [I0l[22]. 

The strong size dependence observed in Fig. [6] for SF 
networks (7 = 2.25 and 2.75) but not for the homoge- 
neous one (7 = 3.25) calls for an anomalous dependence 
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A(N/g) 



FIG. 7: (Color on-line) Collapses of the nume rica l ME solu- 
tion using the anomalous scaling function Eq. (13 1. Densities 
obtained for network sizes iV = 16 x 10**, 64 x 10*, 256 x lO", 
and 1024 x 10* are shown. Data were shifted to avoid overlaps. 



of p on g, as pointed out in the Langevin 10 and random 
walk mapping [22| approaches. Our numerical approach 
allows a more detailed investigation of this issue. Inspired 
by the anomalous FSS prediction [TO] . 



p ^ A/.g, A > ^i/N, 



(28) 



we have performed a further test, analyzing the behavior 
of the QS density, rescaled as g'^p, as a function of N, as 
shown in Fig. |8] The theoretical exponent 6=1, which 
corresponds to the anomalous scaling in Eq. (28), does 
not show up as a plateau in the plots of g''p versus N. 
Interestingly, plateaus are observed if an exponent 6 < 1 
is used instead. For A = 1.040, plateaus are observed for 
b = 0.88 and h = 0.53 for 7 = 2.25 and 7 = 2.75, respec- 
tively. In an analysis performed closer to critical point, 
for A = 1.004, the plateaus are observed with larger ex- 
ponents b = 0.95 and b = 0.56 for 7 = 2.25 and 7 — 2.75, 
respectively. Notice that a scaling consistent with b — 1 
was obtained for 7 = 2.25 but not for 7 = 2.75. Actually, 
the anomalous scaling ([28]), derived from Eq. ( 13 1, is valid 
for ^Jg/N < A < g{k)/kc. For w = 2, the right and left 
sides of this inequality scale as g{k)/kc ~ c^A''"^'^"^)/^ 
and \/gjN ~ c^N~^''^'^'^^^ , respectively, where c-y — 



fcJ(7-2)/(3-7) and( 
are constants of the same order and 



l)(7-2)V(3-7) 
> c^. If 7 is 



close to 3, the exponents involved in the lower and upper 
bounds of A are very close and we cannot make A suf- 
ficiently small to fulfill the upper bound and still larger 
than the lower one, except for numerically unaccessible 
large systems. Therefore, this anomalous scaling can be 
clearly seen only for 7 close to 2. The scaling forms with 
6 < 1 are thus metastable crossovers between the regimes 
p ^ (giV)^^/^ and p ~ A/5, that can last for decades, 
and could only be resolved by simulations in much larger 
systems sizes than those considered in this work (up to 
= 10^ in the ME solutions). 




FIG. 8: (Color on-line) Probing the anomalous FSS A/g'' 
in the ME equation iterations of the supercritical CP for A = 
1.040 (top) and 1.004 (bottom) and a cutoff exponent u — 2. 
Dashed horizontal lines are guides to the eye. 



VII. 



CONCLUDING REMARKS 



The contact process (CP) on scale-free networks shows 
remarkably rich features, even in the simple case of ran- 
dom annealed topologies. In the present work we have ex- 
plored the quasi-stationary (QS) properties of this prob- 
lem by combining QS numerical simulations and a master 
equation approach applied to an approximate mapping to 
a one-step process. The resulting master equation, apart 
from providing quite accurate analytical approximations 
for the asymptotic shape of the QS activity distribution 
at the critical point, can be very efficiently solved numer- 
ically. The QS distribution and the relevant QS quanti- 
ties (density of active sites and characteristic time) de- 
termined in this way show an excellent agreement with 
direct QS numerical simulations of the contact process, 
both at criticality and in the supercritical regime. 

The high accuracy of our data allows to identify strong 
corrections to the scaling in the critical quantities that 
mask the correct finite size scaling exponents obtained 
analytically by means of an exact mean-field solution. 
Both critical density and characteristic time show tenu- 
ous curvatures as functions of the network size N due to 
finite size corrections to scaling that may provide incor- 
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rect exponents if a simple power law decay is assumed. In 
annealed networks, for which the critical point is exactly 
known, we can determine the corrections to scaling ana- 
lytically and thus recover the theoretical exponents in the 
finite size analysis, including the abrupt change when the 
network loses its SF property. The analysis of the super- 
critical region, on the other hand, hints that those finite 
size corrections are also relevant for very large network 
sizes. Indeed, the asymptotic scaling is observable only 
for extremely large values of much larger than those 
possibly attainable with present day computers. 

It is worth noticing that the QS analysis presented in 
this work is equivalent to the Langevin approach devel- 
oped in Ref. in the limit of very low densities (i.e. in 
the critical region, see Appendix) . The main difference is 
that the former starts from an approximation of the orig- 
inal dynamical processes [Eq. ([6])], while the latter repre- 
sents an exact Langevin approach in the coarse-grained 
limit (density approximated by a continuous variable). 
The critical advantage of the present approach lies on 
the fact that our ME analysis allows to determine the 
critical properties in a very intuitive way, as well as to 
easily obtain highly accurate results, free from statistical 
errors, for all quantities of interest. 

Our work opens the path to a more detailed character- 
ization of absorbing phase transitions on scale-free net- 
works in general, and the CP in particular. In the more 
realistic framework of quenched networks, in which edges 
are frozen and do not change, this goal may be however 
hindered by the interplay between corrections to scaling 
and the usual lack of knowledge about the true position of 
the critical point. In fact, the standard characterization 
of the QS state by the usual procedure assuming a sim- 
ple power law of the system size at the critical point may 
be affected by two sources of errors: the analysis may 
be misleadingly done off the critical point and/or be af- 
fected by important scaling corrections. Further work in 
this direction, following the proposed lines, might thus 
help to throw light on the numerical assessment of the 
correct critical scaling of absorbing phase transitions on 
heterogeneous networks. 
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Appendix 

The connection between the Langevin approach de- 
veloped in Ref. [TU] and the one-step process ^ can 
be established by analyzing the respective Fokker-Planck 
(FP) equations. The general form of a FP equation for 
a stochastic variable x is [5D] 

where A{x) and D{x) are the drift and the diffusion 
terms, respectively. The Langevin analysis in Ref. |10| 
yielded A{n) = n[\ - 1 - \e{n/N)] and D{n) = 
2\nA{n/N), where 8 is given by Eq. ([16]) and 

In turn, the drift and diffusion terms of the FP equa- 
tion for an arbitrary one-step process is given by A(n) = 
W{n+1, n)-W{n-l, n) and D{n) = W{n+1, n)+W{n- 
1, n) [30], respectively. Equation (|6| results in exactly the 
same drift obtained in the Langevin approach, while the 
diffusion term takes the form I?(rt) = (l-l-A)n-f A8(7i/A^). 
At low densities, one can expand 8 and A to leading order 
and obtain D{n) ~ 2?t. + 0(8) in both cases. Therefore, 
Langevin and ME approaches are equivalent in the low 
density limit. 
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